#include<iostream>
#include<fstream>
#include<string>
extern "C"{
#include<math.h>
}
#include"Parton.h"
using namespace std;

Parton::Parton(unsigned int Z_temp=1, unsigned int A_temp=1)
:Z(Z_temp),A(A_temp)            //Z and A are parameters for a nuclei
{
	char filename[50];
	ifstream datain;
	double x, Q2;
	unsigned int i, j;
	xMax=61;
	Q2Max=32;
	flavorMax=6;
	lnxstep=log(10)/((xMax-1)/6);
	lnQ2step=log(2.0);
	grid = new double[Q2Max*xMax*flavorMax];
	//read grid data for interpolation
	sprintf(filename,"grid_%d_%d.dat",Z,A);
	cout<<"    Using "<<filename<<endl;
	datain.open(filename);
	for(i=0;i<Q2Max;i++)
	{
		for(j=0;j<xMax;j++)
		{
			datain>>Q2>>x>>grid[(xMax*i+j)*6+0]>>grid[(xMax*i+j)*6+1]>>grid[(xMax*i+j)*6+2]>>grid[(xMax*i+j)*6+3]>>grid[(xMax*i+j)*6+4]>>grid[(xMax*i+j)*6+5];
		}
	}
	datain.close();
}

Parton::~Parton(void)
{
	delete[] grid;
}

//linear interpolation method
double Parton::fitQuadratic(double x, double* px, double* pf) const
{
	double f01=(pf[1]-pf[0])/(px[1]-px[0]);
	double f12=(pf[2]-pf[1])/(px[2]-px[1]);
	double f012=(f12-f01)/(px[2]-px[0]);
	return pf[0]+f01*(x-px[0])+f012*(x-px[0])*(x-px[1]);
}

//quadratic interpolation method
double Parton::fitLinear(double x, double* px, double* pf) const
{
	double f01=(pf[1]-pf[0])/(px[1]-px[0]);
	return pf[0]+f01*(x-px[0]);
}
